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Abstract 

Cluster algorithms for classical and quantum spin systems are dis- 
cussed. In particular, the cluster algorithm is applied to classical O(N) 
lattice actions containing interactions of more than two spins. The perfor- 
mance of the multi-cluster and single-cluster methods, and of the standard 
and improved estimators are compared. 
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1 Introduction 



Simulations of statistical systems near to their critical points is usually a very 
difficult problem, because the dynamical evolution of the system in Monte Carlo 
(or computer-) time slows down considerably. The phenomenon is called critical 
slowing down (CSD). It also occurs in real systems (in real time) and is an inter- 
esting subject to study. In Monte Carlo (MC) simulations, however, where one 
is interested in quantities averaged over the equilibrium probability distribution, 
CSD could lead to a tremendous waste of computer time. The evolution of a 
system is characterized by the autocorrelation time — this defines the rate how 
the system looses memory of a previous state, or in other words, the number of 
MC steps needed to generate a statistically independent new configuration. 

Let us consider an observable A (a function of spins if we consider a spin 
system) and denote by A t its value at a given MC time t. We define the 
autocorrelation time through the average over the equilibrium distribution as 

C AA (t) = (A s A s+t ) - (A) 2 (1) 

For large t this decays exponentially, CAA(t) oc exp(— t/T exPt A) where T exPt A-> 
the exponential autocorrelation time corresponds to the slowest mode in the MC 
dynamics. It is also useful to define a slightly different quantity, the integrated 
autocorrelation time 

+00 -. +00 

T int ,A C A A(t)/C A A(0) = 2 + E C AA(t)/C AA (0) (2) 

t— — 00 t—1 

This is the quantity which is relevant for the statistical error in a MC simulation. 
One usually estimates (A) from the average of n subsequent measurements, i.e. 
through the quantity 

n 

A=-J> (3) 

t=i 

Its statistical error, 5A for n 3> Ti nt .A is given by 



5 A = v /2r i „ fiJ 4 • 5A name (4) 

where 

SA naive = J-C A a(0) (5) 
V n 

Equation (^) means that in order to reach a given accuracy one has to spend 
a computer time oc Ti nt . Unlike the equilibrium averages (A) or Caa(0), the 
autocorrelation time depends on the actual MC algorithm. In general, near a 
critical point (i.e. for £ 3> 1) the autocorrelation time diverges as 

r a c£* (6) 
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where z is the dynamical critical exponent] 



For a local update, as the standard Metropolis algorithm, one has z « 2. 
Some 2d systems reveal their critical properties only at quite large correlation 
length £. In two dimensions the memory allows to consider systems with linear 
size L « 1000 and a correlation length £ k, 100 — 200. For a local algorithm, how- 
ever, (^J) predicts t ~ 10000, which makes the simulation with a local Metropolis 
algorithm hopeless. 

The t oc £ 2 behaviour reminds one of the the random walk — a change at 
some site will propagate to a distance of ~ £ due to random local updating steps 
in a time proportional to £ 2 . 

There are several ways to change the MC dynamics in order to decrease 
the value of z. One way is to perform a deterministic local update instead of 
a random one, staying on the energy surface -E=const, i.e. making a micro- 
canonical step. This over-relaxation algorithm Q can lower the value of the 
dynamical critical exponent down to z ~ 1. 

The other class of MC dynamics updates a collective mode instead of a 
single local variable. In some optimal cases one can reach z 0, i.e. completely 
eliminate (at least from a practical point of view) the problem of CSD. Two 
important methods of this type are 

• Multi-grid algorithms 

• Cluster algorithms 



In the multi-grid algorithms one introduces in addition to the original (fine 
grid) lattice a sequence of lattices with lattice spacings a' — 2a, a" = 4a, 
. . . and updates the corresponding block of spins, i.e. regions of given sizes 
and shapes. For the O(N) vector model the multi-grid algorithm reduces the 
dynamical critical exponent to z » 0.5 — 0.7 ||. As we shall see, the cluster 
algorithms - the topics of these notes - work better for this system. One should 
keep in mind, however, that for other systems the presently available cluster 
algorithms are not efficient, while the multi-grid algorithms can be still used. 
In these lecture notes I discuss the cluster algorithms for classical and quantum 
spin systems. 



3 For more details and for references on statistical time-series analysis see e.g. 
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2 Cluster Algorithm for Classical Spin Systems 



The cluster MC method has been suggested by Swendsen and Wang 0] for the 
Ising model, based on an earlier observation that the partition function of the 
system can be written as a sum over cluster distributions ||. In Ref. || the 
algorithm has been generalized to spin systems, using more general considera- 
tions. Here I follow this derivation since it will be straightforward to apply it 
to general O(N) lattice actions. 

Consider a general O(N) vector model where the energy of the configuration 
S = {S„} is given byQ 

E(S) =J2 Ei(S) ■ (7) 
i 

Here I denotes a given link, a pair of sites I = (n\, n-i) and the interaction term 
Ei(S) depends only on the corresponding spins, Ei(S) — Ei(S ni , S„ 2 ). At this 
point let us extend the class of actions considered in Ref. |(| and consider the case 
of 3-spin (4-spin, . . . ) interaction terms. Let I be a (hyper)link I = (n 1; ri2> ^3) 
and the corresponding interaction term Ei(S) = (S m , S„ 2 , S„ 3 ) (and similarly 
for 4-spin, ...interactions). We shall also assume that the interaction has a 
global symmetry, 

^(<?Si,<?S 2 ,...) = #;(Si,S 2 ,...) , (8) 

for g £ G where G is the corresponding symmetry group. To construct the 
clusters we connect the sites belonging to some link I by a 'bond' with some 
probability pi (S) depending on the spins associated with the link. We assume 
it to be globally invariant as well: 

pi(gSi,gS 2) ...) =pi(Si,S 2 ,...) . (9) 
The probability to produce a given configuration of bonds B is given by 

W (B\S) = Y[p l (S)H(l- Pl (S)) . (10) 

zes l^B 

The next step is to build the clusters. The cluster is the set of sites which can 
be visited from each other going through bonds. We propose now the following 
change in the configuration: spins within a given cluster are transformed glob- 
ally, by some gi G G, i — 1, . . . , A*c, where Nc is the actual number of clusters 
defined by the bond configuration. (Note that different bond configurations can 
result in the same cluster configuration. In fact, only the cluster configuration 

4 For simplicity, we include the factor /3 = 1/kT in the definition of E, i.e. the Boltzmann 
factor is exp(— E) 
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matters here.) The equation of detailed balance (assuming the same a priori 
probability for g~ x as for gi) is given by 

e- E ^w(B\S)w(S -» S') = e- E ^ s '^w(B\S')w(S' S) . (11) 

Here w(S — > S") and u>(S" — > S 1 ) are the corresponding acceptance probabilities 
used to correct the equation (see Ref. Jjj). Because of the global symmetry (^,^|), 
the contributions coming from a link I with all its sites belonging to the same 
cluster cancel and we obtain 

]Je- E ^(l-p l (S))w(S^S') = ]Je- E ^(l-p l (S'))w(S' -> S) . (12) 

l£B l£B 

Let us introduce a convenient parametrization for pi(S): 

Pl (S) = l-e B *W-QiW . (13) 

The fact that pi(S) is a probability requires 

Qi(S)>Et(S) . (14) 

Equation ( |l2"| ) gives then 

[] e-^w(S -5 ) = n e- Q '( s '^(5' - 5) . (15) 

Here in fact only those links contribute which would connect different clusters, 
i.e. are on the 'surface' of clusters. We have succeeded to map our original spin 
system onto a system of Nq 'sites' (representing the clusters 1, . . . , Nq) with the 
dynamical variables gi,i = 1, . . . ,Nq- The interaction between the clusters is 
given by Qi(S) where I is a (hyper)link connecting two or more clusters. (Note 
that by taking Qi(S) — Ei(S), i.e. pi(S) = 0, each lattice site becomes a cluster 
by itself and we recover the local MC method.) 

Let us specify the bond probabilities, i.e. Qi(S). We restrict the proposed 
transformations g to a subgroup H C G to be specified later, and let Qi(S) be 
the maximal energy Ei(S') where the elements g £ H are applied independently 
to all possible sites belonging to I, 

Qi(Si,S 2 , ■ ■ ■) = max E(g 1 S 1 ,g 2 S 2 ,...) ■ (16) 

9i,92,— £H 

Clearly, we have Qi{giSi, g 2 S 2 , ■ ■ •) = Qi(Si, S 2 , . . .) for any gi,g 2 , . . . £ H. For 
this choice of the bond probabilities the remaining factors / ^ B in ( ^5|) cancel 
and one obtains 

w(S -> S r ) = w(S' -» S) . (17) 
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This is a remarkable fact: the accept-reject step needed usually to correct the 
equation of detailed balance is unnecessary, the clusters can be updated inde- 
pendently. This is possible because the clusters are built dynamically i.e. they 
are sensitive to the interaction and to the actual configuration. With a fixed a 
priori given shape of clusters one cannot achieve this. 

The construction is still quite general, but there is a hidden subtle point 
here. If Qi — Ei is too large then the bond probability becomes too large as 
well. As a consequence, the largest cluster will usually occupy almost the whole 
lattice, leaving out only a few small isolated clusters. Although the algorithm 
is still correct, it is also useless in this case: applying a global transformation 
to the whole lattice does not change the relative orientation of the spins, so one 
would update effectively a few small clusters with the computation being done 
on the whole lattice. Before proceeding with the general case of the O(N) spin 
model, let us consider the Ising model as a special case. 

E = — J SnSn+p, , (18) 

where J > and S n = ±1. The global symmetry transformation in this case is 
S n — > gS n , g = ±1. Equation (|l6|) gives 



US) = J , (19) 



and consequently 



, o \ _ i 1 - c 2,7 for parallel spins, 5i5 2 = +1 , / 9f) -* 
Pi^i, b 2 ) - i Q for antiparallel spinSj SlS2 = _ : \ M ) 



This is the bond probability for the Swendsen-Wang algorithm. An im- 
portant observation is that the clusters cannot grow too large: since the bond 
probability is zero between antiparallel spins, all the spins in a cluster will have 
the same sign - the clusters cannot grow larger than the region of same-sign 
spins. In the unbroken phase the latter does not exceed half of the total number 
of sites, hence the main danger of the cluster algorithms - that one ends up with 
a single cluster occupying nearly the whole volume - is avoided in this case. 

Since the clusters do not interact (c.f. (|l7|)) one can choose any of the 
2 N c s ig n assignments with equal probability. One can average over these 2 Nc 
possibilities without actually doing the updates - by introducing the improved 
estimators. For the correlation function one has 

(S x S y ) =J2p(C)(S*S v )c , (21) 
c 



G 



where C is a given cluster distribution appearing with probability p(C) and 
(S x S y )c is the average over the 2 Nc possibilities for C. The latter is trivial, 

(S x S y )c = S xy (C) , (22) 

where S xy (C) is 1 if the two sites belong to the same cluster, and otherwise. 
Hence we have 

(S a S v ) = (6 xy (C)) . (23) 

The improvement comes from the fact that at large distances \x — y\ the small 
value of correlator (S x S y ) <C 1 is obtained by averaging over values +1 and — 1 
in the standard measurement while in the case of the improved estimator S xy (C) 
it comes from averaging over +1 and 0. In other words, the value of Caa(0) in 
is much smaller for the improved estimator than for the original one since 

{(S x S y ) 2 ) = i, <(MC)) 2 > = (MO> « i ■ (24) 

As a consequence, one expects to gain an additional factor 1/ (S x S y ) in computer 
time. This argument is, unfortunately, not complete. The point is that one 
usually measures the correlation function {S x S y ) from an average over the whole 
lattice with \x — y\ fixed. When the site x runs over the lattice the quantity 
S x S y fluctuates more independently than the improved estimator S xy (C) - the 
values of the latter are strongly correlated when large clusters are present. 

A variant of this algorithm has been introduced by Wolff || - the single- 
cluster algorithm. He suggested to build just one cluster starting from a random 
site, reflect all spins in this cluster and start the procedure again. The difference 
between this and the original Swendsen-Wang multi-cluster algorithm is that 
one enhances the probability of updating large clusters since a cluster of size |C| 
is hit with the probability \C\/V, where V is the total number of lattice sites. 
Large clusters evolve still too slowly compared to small ones in the multi-cluster 
method - the single-cluster version corrects for this by updating more often the 
large clusters. (Obviously, one can try to vary the maximal number of clusters 
to be updated within a given MC run, and optimize in the distribution of this 
number.) 

The improved estimator can be modified for the case of the single-cluster 
algorithm ||. For example, the correlation function is given by 

i v l Nc \c-\ i \ mc / i \ sc 

. <25 > 

where the superscripts mc, sc refer to multi-cluster and single-cluster algorithms, 
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respectively. The susceptibility in the unbroken phase is given by 



^(S a Sy) = /^E E 6 xy (C))= (f^^lQl) = (\C\T . (26) 




y * vi \y 

x,y \ i=l a;,j/gC 

The average size of a randomly chosen cluster diverges together with the sus- 
ceptibility when one approaches the critical point - the MC dynamics adjusts 
itself to the long range correlations present in the equilibrium situation. 

In applying the previous considerations to the O(N) spin model one has to 
choose the proper symmetry transformation H in (|l6|). Choosing this to be the 
rotation group - the choice made in || - leads to a bond probability which 
is too large, the largest cluster tends to occupy the whole lattice. It has been 
suggested there to make only small rotations and consider the case of interacting 



clusters as in ( |15[ ). A much better solution suggested by Wolff g is to choose 
the subgroup of reflections with respect to some given random direction r. (The 
derivation in Ref. || was based on embedding an Ising model into the O(N) 
model and applying to it the Swendsen-Wang cluster algorithm. The advantage 
of the present approach shows up in applications to O(N) lattice actions of 
more general type.) The corresponding subgroup contains the identity and the 
reflection of the parallel component of the spin to the vector r, i.e. 

g : S|j — > — $]! , — > iS^ - . (27) 

For the standard O(N) action we have 

and only the first term is affected by the update. Essentially we obtained an 
embedded Ising model with space-dependent ferromagnetic couplings, 

E = - Jn^n^n+p. > with J n ^ = J|5j|5|| + - 1 , e„ = signS 1 !! . (29) 

The corresponding bond probabilities are given by Qi = Ji- Again, because of 
the ferromagnetic nature of the couplings, the regions with e„ > and e„ < 
are not connected by bonds and the size of the clusters is bounded by the size of 
the corresponding regions. Consequently (at least in the unbroken phase) one 
is safe from having clusters with |C| ~ V. For the O(N) vector model both the 
single-cluster O] and multi-cluster p2] method eliminates practically the CSD. 
For this model the cluster algorithm works even better than for the Ising model 
where some CSD is still observed. 

Note that the present formulation can be easily applied to more general 
O(N) actions. Let us consider first the Symanzik (tree-level) improved O(N) 
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action Q, 



'n+jj, 




'nJn+2ft 



(30) 



Here we have two types of links, I — (n,n + /}) and I — (n,n + 2/t), with 

Qi = fl^Ji^l+zil and Qi = T2\SnSl + 2fi\, respectively. A potential danger is 
related to the fact that the second term in ( p0[ ) is anti-ferromagnetic - the bonds 

— II 

of the second type can connect sites with opposite sign of Sti- Fortunately, the 
coefficient 1/12 is sufficiently small, and the largest cluster does not grow too 



Another example of practical importance is provided by the classically per- 
fect lattice action (or the fixed point action of a renormalization group trans- 
formation) fill which has a form 



This lattice action is used to minimize the lattice artifacts (the discretization 
errors). Although introducing the signs e„ = sign(5 , |[) does not turn this action 
into an Ising model (couplings between more than two spins are also present) 
the discussion presented here readily applies to this general action, and the 
procedure eliminates practically the CSD. 

One can also introduce an external magnetic field. The simplest way to 
apply the cluster algorithm is to consider the external field as an extra spin 
which is coupled to all spins, and consider this interaction term on the same 
footing as all other terms in [jl5[| . The cluster to which the external spin 
belongs does not need to be updated, otherwise the procedure is unchanged. 

When modifying the algorithm, the condition of detailed balance has to be 
rechecked carefully. To illustrate this, let us consider the 3d 0(3) model in the 
broken phase. If the random direction f points towards the magnetization M, 
the effective couplings ■/|S'nS , J| + n | become too large, consequently also the size 
of the largest cluster. (This is expected physically - flipping half of the spins 
along the total magnetization does not produce a typical configuration.) Could 
one modify the algorithm by restricting r to be orthogonal to Ml It is easy 
to see that this is not correct. The component M (orthogonal to r) remains 
unchanged by the updates while M" changes from zero to some nonzero value 
hence \M\ always grows. In fact, we violated detailed balance by biasing f with 
the direction of the magnetization. (On the other hand, it is allowed to take r 
orthogonal to the direction of a given spin.) 



large G3. 



E(S) = p(r)(l - S n S n+r ) + ^2 c ( r ''i-ri4}(l-S ni S n2 )(l-Sn 3 S n4 ) + ... 




(31) 
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One can also define improved estimators for the O(N) spin model. Since 
only the signs of Si are updated one can use the relation Q 

{Sj y ) = N(Slsl) , (32) 

and define the improved estimator as 

N\Slsl\5 xy (C) . (33) 

Unfortunately, this procedure introduces an unwanted noise in the correlator 
when S x and S y are correlated significantly (i.e. for \x — y\ < £). To illustrate 
this let us consider S X S X : this is, of course, exactly 1 when measured directly, it 
has no statistical error while the estimator NS x sl fluctuates. The problem is 
cured easily JljJ : one should only choose N random directions fi , r% , ■ ■ ■ , fiv in 
a sequence forming an orthogonal system and apply the cluster algorithm using 
r = fi for i = 1, . . . ,N, otherwise everything remains unchanged, including the 
improved estimator. Unfortunately, for the single-cluster update this trick is 
not very convenient. To cure the problem with the improved measurement and 
keep the advantage that the single-cluster method updates more efficiently the 
large clusters one can do the following: 

1. make n s single-cluster steps, but do not measure anything, 

2. make N consecutive multi-cluster steps with a random basis and perform 
a measurement after each step. 

One can try to optimize the value of n s to get a uniformly small autocorrelation 
time at all scales or even try to vary the maximal number of clusters (1,2, . . .) 
to be built and updated. We shall refer to this as the 'hybrid' method. To illus- 
trate the performance of the cluster method and of the improved estimators I 
made some runs^ on a 256 2 lattice at (3 = 1.8 which corresponds to a correlation 
length £ = 64.78(15) @. For the hybrid algorithm I took n s = 25 w V/{\C\), 
and measured the correlation function C(x) = (SqS x ) using the standard es- 
timator (averaged over the lattice volume) and the improved estimator. The 
measurements in the single-cluster method were done after each n s updates. 
Table 1 shows the results for C{x) at x = 1, 20, 60, 120, and for the susceptibil- 
ity. To compare the performance the squared errors are shown, normalized to 
20000 sweeps for each case. 

Comparing the standard estimators one sees that the single-cluster method 
updates all scales better, except for the shortest distances. The error squared for 
the hybrid method is somewhat larger because we made there ~ 1 single-cluster 

5 Most of these runs were made after the school, while preparing this written version. 
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method C(l) C(20) C(60) C(120) 



H, imp. 0.68796(3) 0.1973(2) 0.0731(3) 0.0354(4) 3560(20) 



S, st. 


2 ■ 


10" 9 


5 


io- 8 


1 


io- 7 


2 ■ 


io- 7 


6 ■ 


IO 2 


M, st. 


1.6 


• io- 9 


3 


io- 7 


9 


io- 7 


1 • 


IO" 6 


2 • 


10 3 


H, st. 


1.3 


• io- 9 


9 


io- 8 


3 


io- 7 


4- 


io- 7 


7- 


IO 2 


S, imp. 


2 • 


10~ 5 


5 


10~ 6 


2 


10~ 6 


7- 


io- 7 


3 • 


10 3 


M, imp. 


1 • 


io- 9 


1 


io- 7 


3 


io- 7 


4- 


io- 7 


1 • 


10 3 


H, imp. 


1 • 


1Q -9 


5 


io- 8 


1 


io- 7 


1.4 


• io- 7 


4- 


IO 2 


H, imp.* 


7- 


io- 10 


2 


io- 8 


5 


io- 8 


6- 


io- 8 


1.5 


• IO 2 



Table 1: Spin-spin correlation function and susceptibility with different variants 
of the cluster algorithm (single-cluster, multi-cluster and 'hybrid'), measured 
with standard and improved estimators. The first line shows the data, the fol- 
lowing lines the error squared, normalized to 20ksweeps. The last line shows the 
(impractical) choice when the multi-cluster steps are used only in the improved 
estimator but not in updating the configuration. 

sweep (i.e. n s = 25 single-cluster updates) followed by 3 multi-cluster updates 
with measurements. A striking observation is that for the single-cluster method 
the standard estimator produces a smaller error than the improved estimator 
- this is naturally connected with the noise discussed above. For the multi- 
cluster or hybrid method the corresponding improved estimator (which cancels 
this noise by taking a random orthogonal basis) has obviously smaller errors 
than the standard estimator, but the gain is about a factor of 2 — 3, much 
smaller than predicted by the naive argument using the relations in (|24|), At 
x = 60 ~ £ the hybrid version gives an improvement in computer time of a 
factor 3 compared to the multi-cluster version and a factor of 20 compared to 
the single-cluster version with (noisy) improved estimator. It seems that the 
best strategy is to use the single-cluster method with the standard estimator or 
the hybrid version with the improved estimator. Of course, at this correlation 
length all these versions of the cluster method are incomparably better than the 
local Metropolis update. 

To separate the effect of multi-cluster updates from that of the multi-cluster 
improved estimator I made a run where after n s single cluster updates in 3 
consecutive steps (in all 3 orthogonal directions) the multi-clusters were con- 
structed and the corresponding improved estimators were measured but the 
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resulting multi-cluster configurations were not updated. Accordingly, only the 
single-cluster steps were counted in the total number of sweeps. These results 
are shown in the last line, as H,imp.*. Note that this is not a practical choice 
since it takes no time to update the clusters once they were constructed. Al- 
though we did not count the construction of multi-clusters (since it does not 
affects the configurations) it counts, of course, in computer time. 

The cluster method has still many interesting applications for the 0(N) 
vector model [18-23]. As an example, let me mention here the determination of 
the interface tension in the Ising model by Hasenbusch j2^]. He considered the 
Ising model on a 2d strip with different (periodic and anti-periodic) boundary 
conditions in the time direction. The boundary condition has been considered 
as a dynamical variable^, e — 1 for periodic and e = — 1 for anti-periodic b.c. 
The idea in Ref. was that when a line of deleted bondsQ cuts through the 
strip in spatial direction one can flip (with probability 1) the region between 
this line and the boundary (at t = Q) flipping the sign of e at the same time. 
The relative number of cases with e — 1 and — 1 is related to the free energy of 
a kink, i.e. to the interface tension. 

Unfortunately, the cluster algorithm does not work well for other models ]25| , 
at least in the most convenient version with independent clusters. The reason 
is always the same: the largest cluster tends to occupy the whole volume. To 
illustrate this, consider the Ising spin-glass, 

E(S) =-J2 J » S i S J > ( 34 ) 

where the random couplings can be positive or negative. In this case (when there 
are frustrated links) the clusters can grow through the boundaries separating the 
regions with 5 = 1 and S = — 1 , and for the couplings of interest one big cluster 
is formed. One can still give up the requirement of independent clusters and go 
back to (|l5|). By decreasing Qi(S) one can make the clusters smaller - the price 
paid for this is that they start to interact and the acceptance probability may 
be extremely small for clusters of reasonable size. This possibility, however, has 
not been properly investigated yet. 



6 A similar trick has been used in [|24j for SU(3) gauge theory to determine the string 
tension. 

7 Sometimes it is more convenient to say that one puts originally bonds everywhere and 
then deletes them with probability pi = 1 — pi . 
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3 Quantum Spin Systems 



We shall consider the example of the 2d anti- ferromagnetic quantum Heisenberg 
model given by the Hamiltonian 

(35) 



where J > and S x = \o x is a spin operator (a x are Pauli matrices) at site x. 
This old model seems to describe the dynamics of the electron spins within the 
copper-oxygen planes of the La2Cu04 material. (Note that the first high-T c 
super-conductor discovered was La 2 _ a: Ba ;c Cu04 with the doping x w 0.15 pq].) 

To simulate a quantum spin system is much more difficult than a classical 
one because - as we shall see - the MC updates have to satisfy some constraints. 
Based on the work of Evertz, Lana and Marcu p7| w ho developed a loop cluster 
algorithm for vertex models, Wiese and Ying |28| worked out an analogous 
procedure for the Hamiltonian (|35|). I shall follow their derivation, with a few 
insignificant simplifications. In particular, I will consider the Id case, and the 
multi-cluster version instead of the single-cluster one. 

First we decompose the Hamiltonian into non-interacting sub-lattices, H = 
Hi + H 2 , where 

H\ = J S2 n S"2n+l , H% = J Sjn-lSzn > (36) 
n n 

and use the Suzuki-Trotter formula for the partition function: 

Z = Tre-e H = lim Tr \ e -^ e -^ H '] N = lim Tr ( e -^ e -^ ffa • • ■) , 

(37) 

where e = 1/N determines the lattice spacing in the Euclidean time direction. 
By inserting a complete set of eigenstates |+) and |— } of a x between the factors 
exp(— e(3Hi) we obtain a classical system with Ising like variables S[x, t) = ±1 at 
each site of the 2d lattice. To see the structure of the corresponding contribution 
it is sufficient to consider a single interaction term in (|36"|), h = J (s a Sf, — j\ 

where the term 1/4 is subtracted for convenience. Since S a Sb — \ — \{S a + 
Sb) 2 — 1, the eigenvalues of h are — J for the singlet state and for the triplet 
state of the total angular momentum S a + Sb- This gives the following transition 
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amplitudes^] for the operator A = exp(—e/3h): 

(+ + \A\ + +) = (--\A\--)=w 1 = l , 

(+-\A\ + -) = (- + \A\-+)=w 2 = Ue^ J + l) , (38) 
(+ - \A\ -+) = (- + \A\ +-)= w 3 = § {e<f» 1) . 

All other matrix elements are zero: 

(+ + \A\ + -} =w 4 = Q (39) 

The separation H = Hi + H 2 leads to a checker-board type interaction of the 
corresponding classical spin system - only the four spin on the corners of 'black' 
squares interact, producing the corresponding factor Wi in the Boltzmann factor. 
The partition function is given by 

ry nAS) n 2 (S) n 3 (S) n 4 (S) / An\ 

Z = y, W l w 2 w 3 w 4 ■■■ (40) 

{S} 

where rii(S) is the number of 'black' plaquettes of type i. Clearly, only those 
configurations contribute for which 714(5) = 715(5) = . . . = 0, i.e. the configu- 
rations should satisfy the corresponding constraints. This requirement causes a 
serious problem for a local updating procedure since most steps are forbidden. 
In the loop cluster approach this problem is avoided - one builds closed loops 
of bonds on the configuration, and flips all spins along a loop. The new config- 
uration automatically satisfies the constraints, and the bond probabilities are 
chosen to satisfy detailed balance. 

The first step is to connect the sites of a 'black' square by bonds with some 
probabilities depending on the type of the configuration. The prescription is 
given in Table 2 (where the time axis is the vertical one): 

It is easy to see that on a finite periodic lattice these bonds form a set of 
closed loops. Flipping all spins along a closed loop leads to a new admissible 
configuration (i.e. no forbidden 'black' plaquettes of type i = 4, 5, . . . appear). 
By flipping the spins on one of the bonds of a plaquette, the transition prob- 
abilities are p(l — ► 2) = 1, p(2 — > 1) = l/w 2 = Wi/w 2 , p(2 — > 3) = W3/1V2, 
p(3 — > 2) = 1. These probabilities satisfy the relations 

w lP (l -f 2) = w 2 p(2 -> 1) , w 2 p(2 -> 3) = w 3P (3 -» 2) , (41) 

and consequently the condition of detailed balance for the equilibrium distribu- 
tion corresponding to (pL0|). This technique made it possible to determine the 

8 A direct calculation leads to a negative number, — U13 for the third amplitude. Redefining 
the sign of the eigenvector |— } at every odd site makes this amplitude positive without affecting 
the others. 
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type 


bonds prob. 


type 


bonds 


prob. 


1 


+ + 

1 

+ + 


2 




p=l/w 2 


3 


-—— + 

1 


2 


+ • 


p' = w 3 /w 2 




+ — — - 




+ • 





Table 2: Bond probabilities for different spin configurations on a plaquette. 
Note that p + p' = 1 for configurations of type 2. 



low energy effective parameters (as the spin wave velocity and spin stiffness) of 
the anti- ferromagnetic quantum Heisenberg model to a high precision |2q] . 

In ( [37|) one has to take the e — > (N — > 00) continuum limit in the time 
direction. Notice that the only non-deterministic choice is for plaquette of type 
2 where the possibility to put the link horizontally (along the spatial direction) 
is 

v ' = ^ 2 = '^Ti = > + 0{ ' 2) ■ (42) 

One can take now a continuum time and reformulate the prescription into a 
continuum language: to turn a path from vertical into horizontal direction in 
time interval dt equals to \ j3Jdt. Modifying the loop cluster algorithm this way 
Beard and Wiese j2i| have shown that this observation allows to design a loop 
cluster algorithm free from discretization error. 

There are still many promising application of the cluster algorithm to quan- 
tum systems but it is beyond our scope to discuss them here. I would like to 
mention only a recent suggestion by Galli M to overcome the 'sign problem'. 



4 Conclusion 



Cluster algorithms work excellently for classical and quantum spin systems. The 
way of choosing the collective modes to be updated depends on the action and 
the actual configuration, and adjusts itself optimally. The Boltzmann weight 
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(in the optimal case) is completely absorbed in the bond probabilities, hence 
the clusters can be updated independently. This leads to an effective update 
of large scale collective modes, to a strong reduction or even elimination of the 
critical slowing down, and to a possibility to introduce improved estimators with 
reduced variance. 
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